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Abstract 

The engineering interest about dune fields is dictated by the their interaction 
with a number of human infrastructures in arid environments. Sand dunes dy¬ 
namics is dictated by wind and its ability to induce sand erosion, transport and 
deposition. A deep understanding of dune aerodynamics serves then to ground 
effective strategies for the protection of human infrastructures from sand, the 
so-called sand mitigation. Because of their simple geometry and their frequent 
occurrence in desert area, transverse sand dunes are usually adopted in literature 
as a benchmark to investigate dune aerodynamics by means of both computa¬ 
tional or experimental approaches, usually in nominally 2D setups. The present 
study aims at evaluating 3D flow features in the wake of a idealised transverse 
dune, if any, under different nominally 2D setup conditions by means of com¬ 
putational simulations and to compare the obtained results with experimental 
measurements available in literature. 

Keywords: dune aerodynamics, Computational Wind Engineering, 3D flow, 
mushroom-like coherent flow structures 


1. Introduction 

The engineering interest about dune fields is dictated by their interaction 
with a number of human infrastructures in arid environments, such as roads and 
railways, pipelines, industrial facilities, farms, buildings (e.g. Alghamdi and Al- 
Kahtani, 2005). Some of such undesired effects are shown in Figure 1. The 
development, shape and migration of dunes depend fundamentally on availabil¬ 
ity of mobile sediment, on the incoming wind directionality, and on the flow 
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Figure 1: Windblown sand interaction with anthropic activities: megadunes surrounding a 
road and farmlands in Dunhuang, Gansu Province, PRC. Photo: I.A. Inman, 2007 (a), linear 
dunes encroaching Nouakchott, capital of Mauritania. Landsat 1565-10032-6, 1974 (b), loose 
sand covering the Aus to Liideritz railway line, Namibia. Photo: K. Dierks, 2003 (c) 


structures of local disturbed wind (’’topographically forced wind” in Geomor¬ 
phology literature), that is on the wind ability to induce sand erosion, transport 
and deposition. In particular, coherent flow structures embedded within fluid 
flow fields are considered to govern the magnitude, form and scaling of sedi¬ 
ment transport events (Bauer et al., 2013), and the dune shape in turn. Hence, 
a wind-focused perspective has been adopted in aeolian dune geomorphology in 
order to investigate the former to explain the latter. In wind engineering, a deep 
understanding of dune aerodynamics serves to ground effective strategies for the 
protection of human infrastructures from sand, the so-called sand mitigation. 
In areas of constant wind direction and under high sand availability, the trans¬ 
verse dune - which has nearly fixed profile in the direction perpendicular to the 
wind - is the prevailing dune type (Livingstone and Warren, 1996). Furthermore, 
ideal sharp-crested transverse dune, i.e. without crest-brink separation (Bauer 
et al., 2013), are usually adopted in fundamental aerodynamic studies for their 
geometric simplicity. Previous studies (e.g. Lancaster, 1995) have revealed that 
transverse dunes generally have upwind (windward, stoss) slope angles ranging 
between 2° < a u < 20° and downwind (leeward) slope angles ranging between 
28° < ad < 34°, i.e. around the sand friction static angle. For sake of clarity, a 
simplified 2D scheme of the transverse dune geometry and of the flow around it 
is given in Figure 2. The aerodynamic behaviour of sand dunes in atmospheric 
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boundary layer belongs to the very general class of bluff bodies and to the gen¬ 
eral one of hills (for a review, see e.g. Bitsuamlak et al., 2004). The incoming 
wind flow is first fairly decelerated at the upwind toe, then strongly accelerated 
along the dune upwind face straight up to crest. The separation of the boundary 
layer then occurs at the crest itself. The naturally inclined downwind face is 
surrounded by a reversed flow region, whose extent is one of the key parameters 
describing the flow. Reattachnrent of the boundary layer occurs far downstream 
the dune crest, being x r the so-called reattachment length. The latter quantity 
is experienced in sharp-edge bluff body aerodynamics to be highly sensitive to 
a number of setup parameters (e.g. Reynolds number, surface roughness, in¬ 
coming turbulence). Much more detailed 2D topological mapping of the wake 
region for isolated transverse dunes under crest-normal flow can be recovered 
in literature (e.g. Walker and Nickling, 2002, 2003; Delgado-Fernandez et ah, 
2011 ). 

While the 2D flow structures above have been extensively scrutinized and re¬ 
viewed in Livingstone et al. (2007), very little is known about three-dimensional 
coherent flow structures in the wake (for an updated review, see Bauer et al., 
2013). On one hand, in the Geomorphology community, such 3D coherent flow 
structures past nominally 2D dunes are usually ascribed to: i. oblique incom¬ 
ing wind, i.e. where a yaw angle ac > 45° gives rise to helical vortices in the 
wake (Allen, 1970; Walker and Nickling, 2002), and/or ii. small crest irregu¬ 
larities, inducing spanwise flow in the wake, spanwise swirling of the reversed 
flow, and secondary recirculation zones (see the recent CFD simulations and 
field measurements in Jackson et al., 2011, 2013; Delgado-Fernandez et al., 
2013, , respectively). It is worth pointing out that from an aerodynamic point 
of view both occurrences involve a 3D setup. On the other hand, to the writers 
best knowledge, 3D coherent structures in the wake of genuine 2D setups (i.e. 
ideal transverse dune under normal wind) remain elusive and scarcely studied 
in literature, if any. Difficulties in both wind tunnel facilities and measurement 
techniques, and in 3D computational simulations are conjectured by the writers 
and other authors (e.g. Walker and Nickling, 2003) to be the cause of such lack 
of knowledge. 

On one side, Wind Tunnel (WT) tests in the aeolian field adopt scaled dune mod¬ 
els spanning across the whole test section width (i.e. S = w), and with rather 
low values of the dune aspect ratio (span/chord ratio S/L), e.g. S/L = 1.28 in 
Walker (2000) and Walker and Nickling (2003), 1.92 < S/L < 10.32 in Dong 
et al. (2007); Qian et al. (2009), S/L = 5.4 in Liu et al. (2011)). Such values of 
the aspect ratios: 

• are by far lower than the ones adopted in other bluff body aerodynamics 
problems characterised by separation, reversed flow and reattachment, 
especially when 3D flow features are expected (see for instance Bruno 
et al., 2014); 

• are suspected to significantly affect the 3D features of the flow in the wake, 
because of the tips effects due to interaction between the boundary layers 
surrounding the dune and the WT side walls. For instance, in Walker and 
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Nickling (2003, Figure 5a, page 1119) the normalised shear stress profile 
along the WT midline (i.e. .46 m far from the WT side wall) does not 
reach nil value in the conjectured reattachment area nor elsewhere in the 
dune wake. 

More generally, difficulties in WT studies are recognised in the proper scaling of 
different setup lengths, i.e. the ones related to the dune geometry, the surface 
roughness, the turbulent length scale of the turbulent incoming flow (Walker 
and Nickling, 2003). 

On the other side, the fundamental computational (CFD) studies on the flow 
field over the transverse dune (e.g. Parsons et al., 2004a, b; Schatz and Herrmann, 
2006; Araujo et al., 2013) usually consider simplified 2D condition (see Living¬ 
stone et al., 2007, for a review). Only recently, Liu et al. (2011) have compared 
the simulated flow around a transverse dune obtained in 2D and ”2.5D” condi¬ 
tions between them and with WT measurements. The ”2.5D” heading indicates 
that despite the computational domain is 3D, its height and width are equal to 
half the WT working section dimensions where reference experimental test are 
performed. In other terms, both left-to-right and bottom-to-top symmetry of 
the flow are conjectured to reduce the domain size and related computational 
costs. Under such an assumption, interesting preliminary results are obtained. 
In particular, 3D flow features of the flow in the near wake of the 2D dune are 
qualitatively observed (Liu et al., 2011, Figure 8a, page 884). According to the 
writers, some questions immediately follow. Which aerodynamic phenomena 
underly and/or trigger the 3D flow features observed in the cited experiments 
and computational simulations? In particular, do the WT side wall effects play 
any role in generating such 3D structures? That is, would similar structures 
having corresponding characteristic lengths obtain if the same dune had, in the 
limit, an infinitely long span (i.e., an infinite aspect ratio)? 

The present study aims at shedding some light on such issues. The sensitiv¬ 
ity of the 3D features of the wake to the setup geometrical scaling is studied. In 
particular, the domain size and the related inlet and side boundary conditions 
are the retained parameters of the study. In such a perspective, the present 
study takes the baton relay from the stimulating study of (Liu et al., 2011) 
conceived in the Geomorphology community, and develops it according to the 
knowledge background of Bluff Body Aerodynamics. The adopted computa¬ 
tional approach is expected to efficiently complement the wind tunnel studies in 
exploring a huge number of setup conditions, where a single parameter is varied 
at the time and border, or even unphysical, conditions can be scrutinized. The 
experimental setup of Liu et al. (2011) is adopted as the reference one. The 
obtained computational results are compared with the experimental ones avail¬ 
able in literature. Finally, a deeper insight in the 3D emerging coherent flow 
structures in the wake is provided. 


4 


2. Wind flow modelling and computational approach 


The incompressible, turbulent, separated, unsteady flow around the dune 
profile is modeled by the classical Time-dependent Reynolds Averaged Navier- 
Stokes (T-RANS) equations, which, in Cartesian coordinates, read: 
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where zlf is the averaged velocity, v! the velocity fluctuating component, p the 
averaged pressure, p the air density and v the air kinematic viscosity. The SST 
k — ui turbulence model first proposed by Menter (1994) and further modified 
in Menter et al. (2003) is used to close the T-RANS equations: 
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where k is the turbulent kinetic energy, ui its specific dissipation rate and u t 
the so-called turbulent kinematic viscosity. The kinetic energy production term 
P k is modeled by introducing a production limiter to prevent the build-up of 
turbulence in stagnation regions: 

~ du ■ 

P k = min (P fc , 10/3* kui) being P k ss 

For sake of conciseness, the definition of the blending function i 7 ) and the values 
of the model constants are omitted herein. Interested readers can find them in 
Menter et al. (2003). The SST k — ui turbulence model is selected for the 
current application because of its proven accuracy in bluff body aerodynamics 
in general (Menter et al., 2003) and in dune aerodynamics in particular (Liu 
et al., 2011). 

At the ground and dune surfaces the so-called sand-grain roughness wall func¬ 
tions are selected for the current application because of their wide use in envi¬ 
ronmental CWE in general (e.g. Blocken et al., 2007a) and the proofs of ade¬ 
quacy obtained in previous 3D simulations of sand dune aerodynamics by Liu 
et al. (2011); Jackson et al. (2011, 2013). In particular, standard wall functions 
(Launder and Spalding, 1974) with roughness modification (Cebeci and Brad¬ 
shaw, 1977) are applied. The required equivalent sand grain roughness height 
is determined as K s = 9.793zo/C s , where C s = 0.5 is the roughness constant, 
and the first cell at wall is limited by n p > Ks, being n p half of the cell size, 
following Blocken et al. (2007b). 

The adopted computational domains and the conditions imposed at their 
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boundaries are the object of a parametrical study. They are detailed in the 
next Section. 

The OpenFoam©Finite Volume open source code is used in the following to 
numerically evaluate the flow-field. The cell-centre values of the variables are 
interpolated at face locations using the second-order Central Difference Scheme 
for the diffusive terms. The convection terms are discretised by means of the 
so-called Limited Linear scheme, a 2nd order accurate bounded Total Varia¬ 
tional Diminishing (TVD) scheme resulting from the application of the Sweby 
limiter (Sweby, 1984) to the central differencing in order to enforce a mono¬ 
tonicity criterion. The pressure-velocity coupling is achieved by means of the 
pressure-implicit PISO algorithm, using a predictor-corrector approach for the 
time discretisation of the momentum equation, whilst enforcing the continuity 
equation. The space discretization is accomplished by a predominantly struc¬ 
tured grid of hexahedral control volumes. Unstructured patterns locally occur 
at the intersection between the surface-fitted grid boundary layer and the carte¬ 
sian grid in the higher part of the domain. The denser grid is located close to 
the wind tunnel floor and side walls, and to the dune surface. The height of the 
control volume adjacent to the wall n w is driven by the sand-grain roughness 
wall function requirements (Blocken et al., 2007b). In particular, n w is set equal 
to n w = 0.2 h in order to satisfy at best both a sufficiently high mesh resolution 
in the normal direction n to the surface, and the requirement n p = n w /2 > Kg, 
being for the adopted mesh n p /Kg = 1.28. The advancement in time is 
accomplished by the implicit two-step second order Backward Differentiation 
Formulae (BDF) method. The adopted time step is equal to At = 0.00025 [s], 
i.e. A tu re f/h = 0.1 dimensionless time unit. 

3. Application setup 

Three setups are adopted and the resulting flow are compared. All of them 
adopt a sharp-crested transverse dune, and generally refer to the experimen¬ 
tal setup described by Liu et al. (2011). Besides such main reference, critical 
comparison is also made to other measurements provided in Dong et al. (2007); 
Qian et al. (2009); Walker and Nickling (2003), the former two being obtained 
for several incoming speeds, the latter in slightly different experimental condi¬ 
tions. Table 1 summarizes the main features of the cited experimental setups, 
where zq is the dune and floor aerodynamic roughness and u re f is the incom¬ 
ing reference speed. It follows that the reference Reynolds number is equal to 
Reh = u re fh/v = 1.7e + 4. 

The setups differ in both the size of the analytical domain in space and in the 
corresponding applied boundary conditions (b.c.). The domain size and the 
kind of b.c. adopted in the setups are depicted in Figure 3, while the profiles of 
the incoming wind (mean velocity u x , turbulence intensity It and length scale 
Lt) are plotted in Figure 4. The setups are briefly commented in the following: 

si only 1/4 of the wind tunnel working section volume is retained. Symmetry 
conditions are imposed on the vertical plane because of the conjectured 
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Table 1: Experimental setups 


Authors 

h 

[mm] 

Ot u 

[deg] 

otd 

[deg] 

L/h 

H 

S/L 

H 

U r ef 

[m/s] 

zo 

[mm] 

Liu et al. (2011) 

Dong et al. (2007) 

25 

10 

30 

7.4 

5.40 

10 

8,10,12,14 

0.1 

Qian et al. (2009) 

Walker and Nickling (2003) 

80 

8 

30 

9.0 

1.28 

8,10,12,14 

8,13,18 

n.a. 


symmetry of the flow. Free stream b.c.s are set on the upper horizontal 
plane to disregard the influence of the WT upper wall. The incoming 
mean wind velocity and turbulence intensity profiles are fitted on the WT 
measurements (Liu et al., 2011, , estimated shear velocity u* = 0.512 
m/s), while turbulence length scale is conjectured constant and equal to 
Lt = 5 mm because of the lack of experimental data (Figure 4). The 
spatial grid involves about 4.e + 5 control volumes. This setup is often 
retained in CFD practice to reduce the number of grid volumes and related 
computational costs when the WT conditions are to be emulated (e.g. in 
Liu et ah, 2011). The setup is intended to discuss the accuracy of such an 
approach; 

s2 the setup exactly reproduces the size of the WT working section in Dong 
et al. (2007); Qian et al. (2009); Liu et al. (2011). No-slip b.c.s are set 
at the four alongwind boundary planes to model the WT walls. A com¬ 
plementary simulation (s2-a, Figure 3) is preliminary performed along an 
empty channel to replicate the WT approaching section described in Liu 
et al. (2011). Periodic conditions at the inlet and outlet of s2-a are set to 
obtain a fully self-developed, horizontally homogeneous incoming bound¬ 
ary layer flow (Blocken et al., 2007b). The resulting u x , It and Lt profiles 
are then imposed at the inlet of the main domain (Figure 4). The spatial 
grid in s2-b involves about 1.5e + 6 control volumes. The setup aims 
at removing the ansatz introduced in si in order to suggest best practice 
guidelines in CFD simulations of WT tests; 

s3 The side walls of the WT working section are replaced by spanwise periodic 
conditions. Two domain spanwise lengths are considered: w^ a = wi is set 
in setup s3a, and w^ a = &W 2 in setup s3b. The height of the domain 
is set equal to 12 h, i.e. larger than the range suggested in Franke et al. 
(2007) 4 h < h\vr < 10 h for external flows in Wind Engineering. The 
spatial grid in s3b involves about 2.7e + 6 control volumes. The setup 
aims at comparing the results to the ones obtained in the s2 setup and 
at evaluating 3D flow features, if any, under 2D, ’’external” incoming 
flow conditions. The profiles of the incoming turbulent length scale and 
turbulence intensity are set in accordance to Richards and Norris (2011) 
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to replicate an external flow. 




no slip | periodic symmetry] | inlet outflow 

Figure 3: Scheme of the adopted setup conditions (not in scale) including b.c. (wind from left 
to right) 


In all setups, Neumann conditions (’’outflow” in Figure 3) involving the velocity 
field and the pressure (null normal component of the stress tensor) as well as k 
and to are imposed at the outlet boundary. No-slip conditions are imposed at the 
floor surface and at the WT walls, if any. Uniform initial conditions are imposed 
at the beginning of the time-dependent simulations. The simulated time for each 
setup is equal to about Tu re f/h = 4000 dimensionless time units, long enough 
to guarantee the convergence of the first and second statistical moments in time 
for all the flow variable, following the convergence check proposed by Bruno 
et al. (2010). The simulations have been performed thanks to the Optiflow 
pc cluster (32 nodes with 4 cores each, 4 Gb RAM per node, Intel Nehalem 
2.8 GHz clock, 1.4 Tflops). An overall cpu time of about 25 hours on 63 cores 
results for the complete parametrical study. 

4. Results 

4-.1. Setup effects on the overall flow regime 

In this section, the time dependent behavior of the simulated flow in the 
three setups is discussed with reference to both bulk and local quantities. The 



































Figure 4: Vertical profiles of the inlet boundary conditions 


former is defined as the aerodynamic coefficient Cpit) of drag force per span- 
wise unit length that results from integration of the stress field on the dune 
downwind face, Figure 5(a). The Power Spectral Density of the drag coefficient 
is plotted versus the dimensionless frequency /* = fh/u re f in Figure 5(b). 
The PSDs in both si and s2 setups clearly show a well defined peak at about 



(a) 



Figure 5: Time dependent flow: PSD of the drag coefficient 


f* ss 0.05 = St, that conversely disappears in s3a and s3b setups (both notated 
”s3” in the following, if not specified otherwise, for the sake of simplicity). At 
such Strouhal number St the power density in setup si is one order of magni¬ 
tude higher than the one in s2, that is in turn by far greater than in s3. The 
high frequency fluctuations (/* > 3e — 1) are due to spurious numerical oscil¬ 
lations of the solution. In summary, the drag force resulting from s3 is almost 
steady ( stdiCo ) ~ 5.6e— 7), while very weak to moderate time fluctuations in 
the remaining setups ( std(C'u) « 7.3e — 5 and std((7u) « l.e — 3 in s2 and 
si, respectively ) suggest a localized unsteadiness of the flow. To prove such 
conjecture, the time history of the longitudinal component of the velocity u x 
is plotted for the three setups and at two probes in the wake (Figure 6). In 
si and s2, the closer the probe to the WT side wall, the larger the velocity 
fluctuations, while in s3 the flow in the wake is almost steady everywhere. An 
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Figure 6: Time dependent flow: spanwise position of the time dependent phenomena 


instantaneous field of the 2 -vorticity on the plane z/h = 1 for s2 is provided in 
the same Figure to shed some light on the underlying fluid flow phenomenon. 
In fact, the interaction between the separated boundary layer at the dune crest 
and the attached one along the WT sidewall induces a 3D flapping along the 
latter. 

As previously introduced, the reattachment length x r is one of the most relevant 
quantities in transverse dune aerodynamics. In the adopted computational ap¬ 
proach, the reattachment point is evaluated as the one in which the longitudinal 
component t x of the shear stress change sign downwind the dune downwind toe. 
Heaving in ming the weak variability in time of the overall flow, the time aver¬ 
aged reattachment point (t — avg{x r )) is evaluated, while its spanwise trend is 
described by the reattachment line, i.e. t — avg(x r (y)). The latter ones obtained 
in the three setups are plotted in Figure 7. 

The following remarks follow: 

• very significant WT wall effects take place close to the dune tips in both 
si and s2 setups. In particular: 

- in si such effects propagate along all the dune span, so that a mean¬ 
ingful estimate of the reattachment length can be obtained by span- 
wise averaging nor by adopting a point wise value, e.g. at the mid 
plane; 

- in s2 the spanwise variability due to the WT side wall is limited 
to a distance from the wall of about d, ss 1.5 L, clearly related to 
the unsteady plumes pointed out in Figure 6. A central segment 
with nearly constant x r can be recognized: at the vertical mid plan 
(y = 0) x r ~ 6.35 h, pretty close to the value measured by Dong et al. 
(2007) (x r ~ 6 h) in the same setup conditions; 

• unexpectedly, a quasi-periodic, ’’festoon-shaped” spanwise reattachment 
line emerges from both simulations s3a and s3b. Such a trend allows to 
estimate some spanwise statistics. 
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Figure 7: Spanwise dependent flow: variation of the reattachment line 


- The first two statistical moments of the reattachment length, i.e. its 
mean value y — avg(t — civg(x r )) and standard deviation y — std(t — 
avg(x r )), are: 

s3a y — civg(t — avg{x r )) = 5.04 h, y — std{t — avg(x r )) = 0.82 h; 
s3b y — avg(t — civg(x r )) = 4.64 h, y — std(t — avg(x r )) = 0.88 h. 

It is worth pointing out that the in both simulations the mean value 
y — avg(t — avg(x r )) is significantly lower than the one obtained at 
the mid-plan in s2. 

- Let us define a characteristic spanwise length (or wavelength) A of 
the reattachment line (see closeup view in Figure 7). Its spanwise 
statistics are: 

s3a y — avg(t — avg{ A)) = 5.61 h, y — std(t — avg(X)) = 0.43/i; 
s3b y — avg(t — avg(X)) = 5.75ft., y — std(t — avg(X)) = 0.36 h. 

The value of the standard deviation y — std(t — avg(X)) is small in 
both simulations, that is the festoon shape is quite stable spanwise. 

- Finally it is worth pointing out that the spanwise dimension w^ a in 
setup s3a accommodates 7 wavelengths (including the two halves at 
the domain side), while the one w S 3 b = 4u> S 3 a in setup s3b accommo¬ 
dates 27 waves plus an oscillation at midspan, i.e. at the far section 
from the lateral periodic b.c.s. 
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In summary, the findings obtained in the Section prove significant effects of the 
WT side walls (si and s2) on the overall flow regime, and suggest that they 
inhibit the triggering of an almost periodic spanwise variability of the flow in 
the wake, as observed in external flow conditions (s3). 

J^.Z. Comparison between WT measurements and CWE results 

In this Section, comparisons between time-averaged WT measurements at 
the dune midspan ( y/L = 0) and the present computational results obtained in 
the different setups are provided. Bearing in mind the quasi-periodic spanwise 
reattachment line in s3 (Fig. 7), also spanwise statistics are evaluated for such 
setup. The z-wise profiles of the longitudinal component of the velocity u x have 



line d - x/h = 0 line e - x/h =1.7 line f - x/h = 10 



(a) (b) (c) (d) (e) (f) 


Figure 8: Horizontal velocity profiles: WT measurements (Liu et al., 2011) and present CWE 
results 


been measured in WT tests by Liu et al. (2011) using PIV at different positions 
in the the vertical mid plan (y = 0) of the WT working section. These profiles 
are compared to the computational ones in Figure 8. The scatter with the WT 
data around the upwind face (Figure 8, lines a-b) are probably due, according 
to Liu et al. (2011), to the limitation of the PIV system, which cannot fully 
resolve the high-speed gradient in the near-surface zone and leads to nonzero 
wind speed at the surface. The agreement in the reversed flow region (Figure 
8, lines e-f) is excellent for setup s2, that replicates the WT conditions. The 
computational results from setup si strongly overestimate the speed deficit in 
the near wake, while the y — avg{t — avg[u x )) in s3 is quite different form the 
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centre plane one and from WT data. In the same setup, a significant spanwise 
deviation occurs in the wake (Figure 8, lines e-f). 

The measured a;-wise profiles of the vertical component of the velocity u y at 
the level z ~ h are reported in Qian et al. (2009). PIV technique has been 



Figure 9: Vertical velocity profiles: WT measurements (Qian et al., 2009) and present CWE 
results 


adopted on the WT vertical mid plan (y = 0, i.e. at a distance from the side 
wall d/L = 2.7) at different reference velocities of the incoming wind. The u y 
profiles obtained in the three computational setups are compared among them 
in Figure 9(a), while the profiles measured at midspan are compared in Figure 
9(b) to the computational ones at different distance from the side wall in setup 
s2. Considerable differences between si and s2 can be observed in the far wake 
only, while a dramatic change in the vertical component of the velocity occurs in 
s3 especially in the near wake. At the dune downwind face (0 < x/h < 1.7 h), u y 
changes its sign from positive (upward flow, si and s2) to negative (downward 
flow, s3). In other terms, the 2D clockwise recirculating region which character¬ 
izes the near wake topology in s2 at mid line no longer holds and it is suspected 
to change in a 3D local flow (see the large spanwise standard deviation of u y ). A 
satisfactory agreement is observed between WT and computational approaches, 
if the same setup and distance from the side wall is adopted (blue points and 
continuous line in Fig. 9-b). Conversely, the closer the vertical plan to the 
side wall, the smaller the upward speed at the upwind face and the higher the 
downward speed at the downwind surface. 
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An analogous post processing is proposed for the z-wise profiles of the shear 
stress magnitude |r| at the floor in Figure 10. The shear stress is a quantity of 



Figure 10: Shear stress profiles: Wind Tunnel measurements (Walker and Nickling, 2003) and 
present CWE results 


particular interest in windblown sand dynamics because above a given threshold 
value Tt it induces the sand grain saltation, i.e. the sand bed erosion. Conversely, 
at |r| < r t sedimentation of the flying grains occurs(Bagnold, 1941). For such 
reason, the |r|/r t obtained by computations in the three setups is plotted in 
Figure 10(a). si fails in predicting the sedimentation just downwind the crest 
and the downwind toe. The lower the blockage ratio, the shorter the sedimen¬ 
tation length l s , i.e. the closer to the dune the point from which erosion takes 
place again. Once more, the spanwise variability of the profile along the wake is 
significant in s3, and the sedimentation length varies from about 0.75L to 1.1L. 
Figure 10(b) compares the shear stress magnitude profiles obtained in s2 with 
the one measured by Walker and Nickling (2003) in their WT test by using 
Irwin-type differential pressure sensors. It is worth pointing out that in this 
case the experimental setup differs from the computational one in the dune 
cross section, in the dune aspect ratio, and in the Reynolds number (see Table 
1). On one hand, the smaller angle of the upwind face a u in WT test induces 
the lower shear stress magnitude along it. On the other hand, the |r| profile 
downwind the crest is mainly affected by the dune aspect ratio, i.e. by the 
distance d of the measurement alignment from the WT side wall. In fact, a 
very good agreement is obtained when the measurements at the WT mid plan 
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(i.e. d = 0.5 S/L ss 0.65) are compared to the computational profile at the same 
distance from the side wall. The closer the alignment to the side wall, the higher 
|r| along the wake. In other terms, in this case study the WT wall side effects 
are by far larger than the ones induced by L/h and/or u re f , even if much more 
emphasis is traditionally given to the latter than to the former. 

4-3. 3D Coherent Flow Structures 

This final Section aims at providing a sound phenomenological reading of the 
spanwise variability in simulated in external flow conditions (setup s3). A deeper 
insight in the flow characteristics is then performed, thanks to the amount of in¬ 
formation available from computational simulations. Advanced post-processing 
and flow visualization techniques are employed to this aim. In particular, the 
velocity and shear stress vector fields are visualized on selected planes by us¬ 
ing the so-called Line Integral Convolution (LIC, Cabral and Leedorn, 1993; 
Stalling and Hege, 1995; Laramee et al., 2003). Some information are provided 
herein, being a technique scarcely used in the Wind Engineering field. 

LIC is based upon locally filtering an input texture along a curved stream line 
segment in a vector field and it is able to depict directional information at high 
spatial resolutions. The directional structure of a vector field can be graphically 
depicted by its stream lines, i.e. paths whose tangent vectors coincide with 
the vector field. Given a pixel Xq of the desired resulting image and numeri¬ 
cally computed a streamline er, passing by Xq = cr(s o), line integral convolution 
consists in calculating the intensity for that pixel as 



(5) 


where s is the arc-length coordinate, T is an input texture, 2 L is the filter length 
and k is a filter kernel normalized to unity. In LIC visualizations of fluid flows, 
a white noise or similar random image is chosen as input texture T. The con¬ 
volution causes pixel intensities to be highly correlated along individual stream 
lines, but independent in directions perpendicular to them. In the resulting 
images the directional structure of the vector field is clearly visible. 

Figure 11(a) and (b) show the LIC visualization of the time-averaged velocity 
field on an ideal inclined plane (/3 = 12°) crossing the near wake from the dune 
crest, and the LIC visualization of the time-averaged shear stress field on the 
dune downwind face and wake floor, respectively. In both figures, the direction 
of the filed is highlighted by superimposed arrowed lines. The most striking 
result is that in Figure 11(a), where the LIC of the time-averaged velocity field 
highlights quasi-steady, ’’mushroom”-like coherent flow structures having a pair 
of owl eyes each. These structures consists of two counterrotating foci with a 
free stagnation point downwind the crest and between the foci. The main di¬ 
mensions of such structures are given in the figure. The corresponding LIC of 
time-averaged shear stress field (Fig. 11-b) is even more complicated and char¬ 
acterized by a ’’skein”-like pattern. The pattern shape appears to be strongly 
related with the mushroom-like structure, especially at the ground floor (i.e. 
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Figure 11: Spanwise 3D coherent structures in the dune wake: time averaged velocity field 
(a), time averaged shear stress field colored by the its modulus (b) 


where the mushroom is closest to the floor). In order to shed some more light 
in such a pattern, it is colored by the modulus of the same shear stress, and 
the spanwise lines along which t v = 0 and t x = 0 (reattachment line, see also 
Fig. 7) are superimposed. For both lines a geometrical reading with respect to 
the LIC pattern can be provided: the former is the locus of points where the 
tangent to the streamlines are aligned with the cc-axis, while the latter is the 
locus of points where the tangent to the streamlines are aligned with the y-axis. 
The mushroom spanwise characteristic length is equal to the wavelength A of 
the ’’festoon” reattachment line. In particular, the spanwise maximum reat¬ 
tachment length correspond to the mushroom center, while the minimum value 
takes place in between two consecutive mushrooms. 

To the authors’ best knowledge, analogous spanwise coherent flow structures 
have been previously recognized at least in two studies involving nominally 2D 
setups. 

Similar spanwise mushroom-shaped flow structures have been visualized by 
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Schewe (2001) around an apparently different aerodynamic setup, i.e. the sep¬ 
arated flow around an airfoil at incidence in the critical regime. In fact, in 
his pioneering and seminal work Schewe has observed and described in detail 
mushroom-like flow structures (Figure 12), combined in different arrangements 
for different values of the airfoil aspect ratio S/L (L = w), and bounded by dis¬ 
turbed regions close to the end plates. On the basis of the experimental results, 



Figure 12: Oil-flow pictures of the mushroom-like flow structures at the upper surface of an 
airfoil, after Schewe (2001) 


Schewe conjectures that: i. the occurrence and number of the mushroom-like 
structures depend on the inclination of the base upper surface; ii. the spanwise 
flow structure remains intact as S/L —> oo; iii. the wavelength A scales with the 
sectional dimension(s) of the airfoil; iv. preferential or more natural spanwise 
arrangements exist, in terms of number of cells and their wavelength. These 
speculations seem to hold, to a certain extent, also in the present case study. 
In fact: i. the downwind face of the dune and the base upper surface of the 
airfoil share the same inclination with respect to the incoming wind (ad ~ 30°); 
ii. setup s3 confirms that the coherent flow structure is qualitatively unchanged 
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by passing from s3a to s3b, where the latter mimics an infinite aspect ratio of 
the dune because of the periodic side b.c. and the increased spanwise length; 
iii. the characteristic length A of the mushrooms found by Schewe (2001) is 
comparable to the one in s3, if both are scaled versus h (A ~ 5.75 h, Fig. 11-a, 
A « 5.6 h, Fig. 12); iv. in both Schewe and present setups, stable arrangements 
are characterized by a given number of spanwise aligned mushrooms. In Schewe, 
an even number of vortices is always observed, but the end-plate disturbances 
considerably extend spanwise; conversely, in s3 an odd number of mushrooms 
emerges in both s3a (7) and s3b (27) setups (see Fig. 7), despite the spanwise 
length in s3b is four time the one in s3a. It is worth pointing out that the 
quasi-periodic mushroom structure is corrupted at midspan, if the dune length 
S = w is not an odd multiple of A (s3b). Despite the analogies above, it is 
worth pointing out that the flow structures in Schewe (2001) are recognized only 
in the critical regime, while the occurrence of the same condition is unlikely in 
the present case study. 

More recently, analogous 3D, global, quasi-steady vortices have been simulated 
by Spalart et al. (2014) in an even more different fluid flow setup, i.e. the nom¬ 
inally 2D, high Reynolds number, fully developed turbulent Couette flow. The 
RANS-based simulations of Spalart et al. (2014) (including the SST k — uj one), 
successfully reproduce the same flow structures previously observed and sim¬ 
ulated in experiments and DNS, respectively. In the present perspective, it is 
worth pointing out that: i. the vortices are arranged in a quasi-periodic z-wise 
system of counter-rotating pairs (Fig. 1, Spalart et al. (2014)), analogously to 
the arrangement in Figure 11(a) ; ii. the distribution of the skin-friction coef¬ 
ficient in the 2 direction (Fig. 6, Spalart et al. (2014)) qualitatively matches 
the ’’festoon-shaped” spanwise reattachment line in Figure 7; iii. the simulated 
flow structures hold in a wide range of Reynolds numbers, and the effect of the 
vortices is only weakly dependent on Re. In their outlook, Spalart et al remark 
that there is a convincing scale separation between the simulated flow structures 
and the turbulence which is represented by the RANS model. The vortices have 
a lateral scale much smaller than the scale of the geometry (which is infinite in 
the setup), and random locations in z, giving them the nature of an instability. 
The remarks above qualitatively hold also in the present case study. Despite 
such analogies, it is worth pointing out that the flow structures in Spalart et al. 
(2014) are somewhat triggered by an initial velocity field mimicking a periodic 
system of counter-rotating vortex pairs in the {y — z) plane, while in the present 
case study the mushroom-like structures spontaneously develop from uniform 
initial conditions. 


5. Conclusions 

The present study points out some emerging 3D coherent flow structures 
in the wake of a transverse dune under different setup conditions by means of 
computational simulations and compares the obtained results with a number of 
experimental wind tunnel measurements available in literature. The comparison 
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shows a good and robust agreement of the CWE results to WT data, and it goes 
beyond: it puts in evidence common misunderstandings in setting up computa¬ 
tional and experimental models for CWE/WT validation. Some good practices 
in dune aerodynamics CWE simulations and WT tests are recommended: i. 
CWE-WT comparisons require the careful simulation of the WT setup, includ¬ 
ing the wind tunnel geometry; ii. forcing the flow symmetry by b.c. in CWE 
can strongly affect the results; iii. high values aspect ratio are recommended in 
WT test to replicate external flow conditions, e.g. w/L 10. 

Surprisingly, emerging mushroom-like, coherent flow structures in the dune wake 
are clearly shown by a deep analysis of the CWE results. They are compared 
to analogous flow structures arising from rather different aerodynamic setups 
belonging to referential literature in fluid dynamics. It follows that the studied 
setup can be ascribed to a wider class of 3D flows occurring under 2D nominal 
conditions. According to the authors such class of flow seems to be of high 
relevance in modern developments of fluid dynamics, and it still remain poorly 
understood in its general features, if any. The authors hope further indepen¬ 
dent studies will appear and further develop their contribution. In particular, 
the precise definition of the aerodynamic regime(s) at which mushrooms grow in 
the dune near wake remain an open issue. An even wider computational study 
would be required to discuss the permanence and/or changes of such coherent 
flow structures form the model scale to the full scale, e.g. by varying Reynolds 
number and/or surface roughness. 
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